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Abstract 

We study an ad hoc extension of the Lattice-Boltzmann method that allows the simulation of 
non-Newtonian fluids described by generalized Newtonian models. We extensively test the accuracy 
of the method for the case of shear-thinning and shear-thickening truncated power-law fluids in 
the parallel plate geometry, and show that the relative error compared to analytical solutions 
decays approximately linear with the lattice resolution. Finally, we also tested the method in 
the reentrant-flow geometry, in which the shear-rate is no-longer a scalar and the presence of two 
singular points requires high accuracy in order to obtain satisfactory resolution in the local stress 
near these points. In this geometry, we also found excellent agreement with the solutions obtained 
by standard finite-element methods, and the agreement improves with higher lattice resolution. 
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I. INTRODUCTION 



Since its origin, more than 15 years ago, the Lattice Boltzmann Method (LBM) has 
proved to be a powerful numerical technique for the simulation of single and multi-phase 
fluid flows in complex geometries. In fact, the LBM has been successfully applied to differ- 
ent problems in fluid dynamics and the interest on the method has grown rapidly in recent 
years. The LBM is particularly suited for complex geometries and interfacial dynamics, 
and its initial applications included transport in porous media, multiphase and multicom- 
ponent fluid flows [J. It was then adapted by Ladd and others to simulate particle-fluid 
suspensions jl2|. It has also been applied to high Reynolds number incompressible flows and 
turbulence, and the implementation of thermal and compressible schemes is being actively 
pursued |3|. One advantage of the LBM is that data communications between nodes is 
always local, which makes the method extremely efficient for large-scale, massively-parallel 
computations (see Ref. J] for an interesting discussion on the LBM capabilities compared 
to the existing continuum-based computational fluid dynamics methods). Another property 
of the LBM that has lately attracted considerable attention is the microscopic origin of 
its mesoscopic kinetic equations, which could therefore readily incorporate molecular level 
interactions. This makes the LBM very compelling for microscale fluid dynamics in microflu- 
idic devices P which typically present non-continuum and surface-dominated effects (e.g. 
high Knudsen number conditions, electrokinetic and wetting phenomena). This microscopic 
based approach also makes the LBM a good candidate for hybrid, multi-scale simulations 
of fluid flows. 

In contrast, the extension of the LBM to non- Newtonian fluids has received limited 
attention so far, in spite of the fact that a reliable extension to the LBM to simulate non- 
Newtonian flows would be very valuable; for instance, in studies of transport in geological 
porous media, an area in which the LBM has been extensively applied jy, E, 0] due to its 
simple implementation in complex geometries, In addition to geological systems, the flow of 
non-Newtonian fluids is commonly found in many areas of science and technology. 

In this work, we study an ad hoc modification of the LBM, first introduced by Aharonov 
and Rothman |9| , in which the local value of the viscosity depends on the strain-rate tensor. 
We show that this modification to the LBM accurately describes the flow of truncated power- 
law fluids, both shear-thinning and shear-thickening, not only in unidirectional flows (parallel 
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plates geometry) but also in two-dimensional flows with simultaneous shear components in 
more than one direction (reentrant corner geometry). 



II. LATTICE-BOLTZMANN METHOD 



The LBM can be viewed as an implementation of the Boltzmann equation on a discrete 
lattice and for a discrete set of velocity distribution functions [lo| . 

/i(x + e,Ax, t + At) = /i(x, t) + n<(/(x, t)), (1) 

where fa is the particle velocity distribution function along the ith direction, f^(/(x, t)) is 
the collision operator which takes into account the rate of change in the distribution function 
due to collisions, Ax and At are the space and time steps discretization, respectively [lj|. 
Then, the density p and momentum density pu are given by the first two moments of the 
distribution functions, 

P = ^2fu pu = ^f i e il (2) 

i i 

where we assumed that the discretization is consistent with the Boltzmann equation, x + e, 
corresponding to the nearest neighbors of the point x. Note that in the previous equation, 
and in the reminder of the article, all quantities are rendered dimensionless using Ax and At 
as the characteristic space and time scales, respectively. Also note that, as we are concerned 
with incompressible flows, we do not need to introduce a dimension of mass. 



A. BGK Approximation 

Assuming that the system is close to equilibrium the collision operator is typically lin- 
earized about a local equilibrium distribution function, f^ q , and assuming further that the 
local particle distribution relaxes to equilibrium with a single characteristic (relaxation) time 
r, we arrive at the Bhatnagar, Gross and Krook (BGK) approximation of the LBM 



a: 



/„(x + e,Ax, t + At) = /,(x, t) + /t(X,t) r(X,t) , (3) 

T 

where the relaxation time r is directly related to the kinematic viscosity of the fluid, v = 
(2r-l)/6. 
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B. Non-Newtonian flows 

The ad hoc extension of the LBM proposed by Aharonov and Rothman P] to simulate 
non-Newtonian fluids consists of determining the value of the relaxation time r locally, 
in such a way that the desired local value of the viscosity is recovered. The viscosity is 
related to the local rate-of-strain through the constitutive equation for the stress tensor. A 
commonly used model of non- Newtonian fluids is the generalized Newtonian model, in which 
the relation between the stress tensor, a^, and the rate-of-strain tensor, Dij, is similar to 
that for Newtonian fluids, er^ = 2fiDij, but with \x a function of the invariants of the local 
rate-of-strain tensor, \x = ujjPij)- I n particular, we are interested in a widely used model: 
the power-law expression fi = m^™" 1 , where n > is a constant characterizing the 
fluid. The case n < 1 correspond to shear-thinning (pseudoplastic) fluids, whereas n > 1 
correspond to shear-thickening (dilatant) fluids, and n = 1 recovers the Newtonian behavior. 
The magnitude of the local shear-rate 7 is related to the second invariant of the rate-of-strain 
tensor, 7 = ^DijDij, where the components of the rate-of-strain tensor, D^, are computed 
locally from the velocity field. In particular, after obtaining the instantaneous velocity field 
from the LBM we then compute from a first-order finite-difference approximation to the 
local derivatives of the velocity. 

However, there is an obvious obstacle to a direct implementation of the power-law fluid 
in the LBM, in that the effective viscosity diverges for zero shear rates (7 = 0) in a shear- 
thinning fluid ( n < 1). Analogously, the viscosity becomes zero for a shear-thickening fluid 
at zero shear rates. In previous studies it is not clear how this problem was avoided. 

Clearly, both limits are unphysical and, in fact, it is known that many non- Newtonian 
fluids exhibit a power-law behavior only in some range of shear-rates, and a constant viscosity 
is observed outside that range Q|. Here, we used the simplest model of such fluids: the 
truncated power-law model, 



u(j) = n(ri)/p 



. (n-l) • , ■ 

m% 7 < 7o 

my {n ^ 7o < 7 < 7oo ( 4 ) 

. (n-l) • . • 

mToo 7 > 7oo 



Using the truncated power-law model has an additional advantage in the LBM. It is 
well known that the LBM can accurately simulate viscous flows only in a limited range of 
kinematic viscosities. The method becomes unstable for relaxation times close to r > 1/2 



12j (small kinematic viscosities, v ^$ 0.001) and its accuracy is very poor for r > 



i Q 



(relatively large kinematic viscosities, v ^ 1/6)- Therefore, we set the lower and upper 
saturation values of the kinematic viscosity in Eq. to v min = 0.001 and v max = 0.1. It is 
clear that the maximum value of the viscosity corresponds to the value at zero shear rate 
for shear-thinning fluids (n < 1), whereas the opposite is true for shear-thickening materials 
(n > 1). Note that setting the value of the maximum model viscosity to v max = 0.1 
for a given maximum fluid viscosity v* max and a spatial resolution Ax simply corresponds 
\o choosing a particular value of the time step in order to satisfy, — v max (Ax 2 /At) 
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. Since the kinematic viscosity scales with Ax 2 /At, to keep the dimensionless viscosity 
constant we shall rescale At according to the previous relationship when we increase the 
number of lattice nodes N, that is, since Ax oc 1/N then At cx 1/N 2 . In what follow we use 
the three-dimensional (face-centered-hypercubic) FCHC-projection model of the LBM with 
19 velocities (D3Q19 following the notation in Ref. 



III. FLOW BETWEEN PARALLEL PLATES 



We first test the proposed LBM for non-Newtonian flows in a simple unidirectional flow, 
the flow between two parallel plates separated a distance b in the ^-direction (Hele-Shaw 
cell) in the presence of a pressure gradient in the x-direction. We use periodic boundary 
conditions in both x and y directions. The resulting flow field is unidirectional, with v x (z) 
the only non-zero velocity component, the rate-of-strain is a scalar function of the local 
velocity, 7 = \dv x /dz\, and the Navier-Stokes equations are greatly simplified jl^ . 

In order to compute the exact solution to the Navier-Stokes equations for a pressure driven 
flow of a truncated power-law fluid in the Hele-Shaw geometry we split the system into (in 
principle) three different regions. We shall describe the regions between z = and z = b/2 
and the analogous regions for z > b/2 follow by symmetry. The first region we consider is 
the high-shear rate region close to the walls, Region if for z < z^, in which the shear rate 
exceeds 7^, and the fluid is Newtonian with effective kinematic viscosity = mjoo , the 
second one is the intermediate region in which the fluid behaves as a power law according 
to Eq. (jlj). Region / for Zh < z < zf, and the last one is the low-shear rate region close to the 
center of the channel, Region L for z\ < z < b/2, in which the shear rate is lower than 70 
and the fluid is again Newtonian, but with kinematic viscosity vq = mj^ 1 . Matching then 
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FIG. 1: Comparison between a Lattice-Boltzmann simulation and the analytical solution for the 
flow between two parallel plates separated a distance b = 10. The power-law exponent of the fluid 
is n = 0.50 (shear-thinning). The pressure gradient is VP = 6 x 10~ 6 , p = 1, vq = 0.1, = 0.001, 
m = 10" 3 , and N = 400. The circles correspond to the Lattice-Boltzmann simulations. The 
solid line corresponds to the analytical solution given by Eq. (fH|l. Also shown, in dashed lines, are 
the continuation of the Newtonian and Power-Law solutions outside their regions of applicability. 
The vertical, dashed lines correspond to the transition points, Z[ and z% = b — Zi, between the low 
shear-rate region L, and the region of intermediate shear rates, /. 



the solution obtained in each region with the conditions of continuity in the velocity and 
the stress, we obtain the general solution to the problem, in terms of the pressure gradient 
Gp = -VP, 



v x (z) 



ra + l 



k (4) z{h ~ z ) + a * 



< z < z h 

Zh< z < Zl 

zi<z< b/2 



(5) 
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with the constants ot\ and «2 given by, 
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Clearly, the number of regions that coexist will depend on the magnitude of the imposed 
pressure gradient G. For very small pressure gradients, G C 1, both z h and z\ become 
negative (see the previous equation), which means that shear rates are smaller than 7 
across the entire gap and only Region L exists. As G increases, there is a range of pressure 
gradients, < (b/2)G < mj^, for which z\ > but Zh < 0, and therefore regions L and 
/coexist. Finally, for G > (2/b)m A f^ we obtain Zh > 0, and all three regions are present in 
the flow. Note that for large values of G both transition points converge to the center of the 
cell, Zi, Zh — > 6/2. Thus, Eq. ((7j) allows us to choose the appropriate value of G in order to 
investigate the different regimes. 

We performed a large number of simulations for different values of the power-law exponent 
n. Specifically, we consider two shear-thinning fluids, n = 0.50 and n = 0.75, and two shear- 
thickening fluids, n = 1.25 and n = 2.00. In all cases we performed simulations for two 
different magnitudes of the external forcing: one for which the region of low shear rates L is 
important, that is relatively small pressure gradients for which z\ ~ b/A; and a second one 
in which the fluid behaves as a power-law fluid almost in the entire gap, that is z\ ~ b/2. In 
both cases the shear-rate does not exceeds 7oo. In Fig.^we present a comparison between the 
Lattice-Boltzmann results and the analytical solution given in Eq. (JHJ) for a shear-thinning 
fluid with power-law exponent n = 0.50. The simulation corresponds to a relatively small 
pressure gradient for which the region of small shear-rates is large, zi ~ 6/4. Both regions, 
Region L in which the fluid behaves as a Newtonian one, and Region I in which the effective 
viscosity is a power-law, are shown. The agreement with the analytical solution is excellent, 
with relative error close to 0.1%. In Fig. El we present a similar comparison between the 
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FIG. 2: Comparison between a Lattice-Boltzmann simulation and the analytical solution for the 
flow between two parallel plates separated a distance b = 10. The power-law exponent of the 
fluid is n = 2.00 (shear-thickening). The pressure gradient is VP = 5 x 10" 6 , p = 1, u = 0.001, 
z^oo = 0.1, m = 10 -3 , and N = 400. The circles correspond to the Lattice-Boltzmann simulations 
and the solid line corresponds to the analytical solution given by Eq. ©. The vertical dashed lines 
correspond to the transition points, z\ and z\ = b — z\, between the low shear-rate region L and 
the intermediate shear rates region I. 

LBM and the analytical solution, but for a shear-thickening fluid (n = 2.00) which behaves 
as a power-law fluid across almost the entire channel. Again the agreement is excellent with 
relative error smaller than 0.1%. 

Finally, for each of these cases we run a series of simulations in which the number of 
lattice nodes, N, in the direction of the gap was increased from 10 to 400, and computed the 
relative error of the LBM results compared to the analytical solutions, defined as ^^(1 — 
v lbm jyAnai.y.^ j n order to obtain the accuracy of the LBM as a function of the number of 
nodes, we simulated the same physical problem but changed Ax from 1 to 0.025. In addition, 
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FIG. 3: Relative error of the LBM compared to the analytical solution for the flow between parallel 
plates, as a function of the number of lattice points used in the simulations. The points correspond 
to simulations with the LBM for four different fluids, two shear-thinning fluids (n = 0.50 and 
n = 0.75) and two shear-thickening fluids (n = 1.25 and n = 2.00). For all fluids we also present 
results corresponding to two different regimes: one at high pressure gradients, in which the low- 
shear-rates region, region L, is small (Power-Law) and the other one at intermediate pressure 
gradients in which both regions L and / are comparable (note that with the exception of n = 0.75 
both regimes give almost exactly the same relative error and, in fact, the corresponding points 
overlap almost entirely). In all cases, we increased the lattice resolution until the relative error was 
on the order of 0.1%. The solid line shows the general trend of the data, 1/N. 

since the accuracy of the LBM depends on the model viscosity, we also changed At according 
to At = Ax 2 , so that the model viscosity remains the same, independent of the number of 
nodes. Then, in order to compare the velocity field always at the same physical time since 
startup, the number of time steps was increased inversely proportional to At (reaching ~ 10 8 
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time steps for N = 400). In Fig. El we present the results obtained for the different fluids 
and different pressure gradients. It is clear that, in all cases, the relative error decreases, 
approximately as 1/N, as the number of nodes is increased, and eventually becomes of the 
order of 0.1% (an arbitrary target accuracy that we set for our simulations). The error was 
found to be independent of the pressure gradient, or the size of the non- Newtonian region 
J, but strongly depends on the power-law exponent. In particular, the relative error seems 
to increase as the magnitude of (1 — n)/n increases, with the error in the Newtonian case 
(n = 1) decaying faster than 1/N. This is probably related to the first-order finite-difference 
approximation used to compute the spatial derivatives of the fluid velocity which determine 
the local viscosity through Eq. (j3J). It would then be possible to improve the accuracy of 
the method by implementing a higher order approximation of the local shear-rates. 

IV. REENTRANT CORNER FLOW 

In the previous section we tested the LBM for non-Newtonian fluids in a Hele-Shaw 
geometry and found excellent agreement with the analytical solutions as the number of 
nodes was increased. In that case, the flow is unidirectional and therefore the shear-rate is 
a scalar, which is a rather simple type of flow. In contrast, we shall now test the LBM in a 
more demanding geometry, that is the reentrant corner geometry sketched in Fig. 0] In this 
case, the shear-rate is no longer a scalar as in the Hele-Shaw geometry and, in addition, the 
presence of two singular points, located at the entrant and reentrant corners (see Fig. HJ), 
requires high accuracy in order to obtain satisfactory stress resolution near these points 
(although no analog to the Moffatt's analysis near the corner is available for non-Newtonian 
fluids, it is believed that not-integrable stress singularities develop in this case, and numerical 
techniques do not always converge ^j]). Motivated by these issues we simulated the flow 
in the reentrant corner geometry using the LBM for a shear-thinning fluid with power-law 
exponent n = 0.50. In Fig. 0] we present the streamlines corresponding to the computed 
velocity field, obtained for a pressure gradient VP = 10 -5 , where the recirculation region 
inside the cavity can be observed (note that the separation between streamlines was chosen 
for visualization purposes only and it is not related to the local flow rate, since the magnitude 
of the fluid velocity sharply decays inside the cavity). The corresponding Reynolds number 
is Re = 4, computed with the maximum viscosity, i>o, and the measured mean velocity. In 
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FIG. 4: Streamlines in the reentrant flow geometry. The flow direction is shown at the center 
of the channel. Note the asymmetry due to inertia effects. The dashed lines show the lines in 
which we compare the solutions of the LBM method with the solutions obtained by finite-element 
calculations. 

fact, the streamlines shown in Fig. 0] are fore-aft asymmetric, due to inertia effects, which 
are absent in low-Reynolds-number flows. We also computed the local magnitude of the 
shear-rate, related to the local stress field through the constitutive relation given by Eq. 
In Fig. El we present a contour plot of the magnitude of the shear-rate in the reentrant 
corner geometry with the lowest level in the contour plot corresponding to 70 . It can be 
seen that the fluid is Newtonian in small regions at the center of the channel and inside the 
recirculation region. It is also clear that, as discussed before, both the entrant and reentrant 
corners are singular points where the shear-rate increases to its highest values in the system. 
Finally, in order to perform a more quantitative test of the LBM, we solved the problem 
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FIG. 5: Contour-plot of the local magnitude of the shear-rate, as computed with the LBM. The 
lowest level of the plot corresponds to 70, that is the fluid behaves as Newtonian in those regions. 
High shear-rates, and accordingly high shear-stresses are localized at the entrant and reentrant 
corners. 

numerically using the finite-element commercial software FIDAP (Fluent Inc.). In Figs. El 
and [7| we compare the solutions obtained with the LBM for different resolutions, ranging 
from N = 40 x 40 to N = 320 x 320, with the finite-element results obtained with FIDAP. 
The comparison is made along two lines, one oriented along the flow direction (dashed line 
at Y = 16 in Fig. HJ) and a second one oriented in the perpendicular direction (dashed line 
at X = 20 in Fig. HJ). m Fig- El we compare the velocity along the channel, u x , in the line 
perpendicular to x that is located at the center of the system (X = 20, see Fig. HJ). A 
velocity profile similar to that in a Hele-Shaw cell is observed for < y < 20, as well as 
the recirculation flow inside the cavity (see the inset). In Fig. [7] we plot the velocity in the 
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FIG. 6: Velocity component along the channel, u x , plotted in a line perpendicular to the flow 
(dashed line X = 20 in Fig. 0J). We compare the results of the finite-element calculations (solid 
line) with the results of the LBM (points) for different lattice resolutions. In the inset we plot the 
velocity profile inside the recirculation region. 

vertical direction, u y , along a horizontal line close to the upper wall of the channel (Y = 16, 
see Fig. HJ). It is clear that there is some fluid penetration into the cavity (entrant flow) 
in the first half of the channel and some reentrant flow in the second half (note that, as 
mentioned before, the flow is not symmetric about x = 20 due to inertia effects). In both 
cases we found an excellent agreement between the two methods, and the agreement clearly 
improved as the number of lattice nodes was increased in the LBM (the number of elements 
in the finite-element computations was fixed to 100 x 100). 
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FIG. 7: Velocity component perpendicular to the flow direction, u y , plotted in a line parallel to 
the flow and close to the top wall of the channel (dashed line Y = 16 in Fig. 0J). We compare 
the results of the finite-element calculations (solid line) with the results of the LBM (points) for 
different lattice resolutions. 



V. CONCLUSIONS 



We have extensively tested an ad hoc modification of the Lattice-Boltzmann Method 
that extends its use to Generalized Newtonian fluids, in which the non-Newtonian character 
of the fluids is modelled as an effective viscosity. Specifically, we calculated the accuracy 
of the method for Truncated Power-Law fluids and showed that the relative error decays 
(linearly) as the resolution of the lattice (number of lattice points) is increased. The error 
was computed directly from the analytical solutions of the problem. The same trend was 
observed for both shear-thinning (n < 1) and shear-thickening (n > 1) fluids, as well as for 
intermediate and high shear-rates. In all cases the relative error was of the order of 0.1% for 
the highest resolution employed. Finally, we also tested the method in the reentrant flow 
geometry and showed that it is in excellent agreement with the solution obtained by means 
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of finite-element calculations. Again, the accuracy of the method was shown to increase 
with the number of lattice points. 
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